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ABSTRACT 

The three Minkowski functionals and the recently defined length of the skeleton are estimated for the 
co-added first-year Wilkinson Microwave Anisotropy Probe ( WMAP) data and compared with 5000 
Monte Carlo simulations, based on Gaussian fluctuations with the a-priori best-fit running-index power 
spectrum and WMAP-like beam and noise properties. Several power spectrum-dependent quantities, 
such as the number of stationary points, the total length of the skeleton, and a spectral parameter, 7, 
are also estimated. While the area and length Minkowski functionals and the length of the skeleton 
show no evidence for departures from the Gaussian hypothesis, the northern hemisphere genus has 
a x 2 that is large at the 95% level for all scales. For the particular smoothing scale of 3?40 FWHM 
it is larger than that found in 99.5% of the simulations. In addition, the WMAP genus for negative 
thresholds in the northern hemisphere has an amplitude that is larger than in the simulations with a 
significance of more than 3 a. On the smallest angular scales considered, the number of extrema in the 
WMAP data is high at the 3 a level. However, this can probably be attributed to the effect of point 
sources. Finally, the spectral parameter 7 is high at the 99% level in the northern Galactic hemisphere, 
while perfectly acceptable in the southern hemisphere. The results provide strong evidence for the 
presence of both non-Gaussian behavior and an unexpected power asymmetry between the northern 
and southern hemispheres in the WMAP data. 

Subject headings: cosmic microwave background - cosmology: observations - methods: statistical 



1. INTRODUCTION 

Over the last few years much interest has been focused 
on constraining non-Gaussian contributions to the cos- 
mic microwave background (CMB). This interest is based 
on three facts: First, if the CMB anisotropy field is in 
fact Gaussian, the power spectrum (or, equivalently, the 
two-point correlation function) contains all the statisti- 
cal information required to characterize the field. Under 
the Gaussian hypothesis, therefore, one may compress 
the statistical information content of the full data set, 
consisting of several millions of data points, into one sin- 
gle vector, Ct, of perhaps a few thousand power spec- 
trum coefficients. Second, most current theories predict 
a Gaussian (or at least nearly Gaussian) anisotropy field, 
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and significant departures from Gaussianity could point 
toward new physics (e.g., Wang & Kamionkowski 2000; 
Lyth & Wands 2002; Bouchet et al. 2002; Bernardeau 
& Uzan 2003). Finally, residual foregrounds or system- 
atic artifacts caused by the instrument or data processing 
arc likely to manifest themselves as non- Gaussian signal 
contributions. Tests for non-Gaussianity afford possible 
methods for detecting, characterizing, and subsequently 
removing such effects. Further impetus for the study of 
such signatures has arisen in recent months from puta- 
tive detections of non-Gaussian features in the Wilkin- 
son Microwave Anisotropy Probe ( WMAP; Bennett et al. 
2003a) first-year sky maps (Coles et al. 2004; Naselsky et 
al. 2003; Vielva et al. 2004; Copi, Huterer, & Starkman 
2004; Park 2004; Eriksen et al. 2004). 

There is no unique way for a random field AT(9, <fi) 
to manifest non-Gaussianity; hence, there is no generi- 
cally superior, high-sensitivity test for all of the possible 
ways in which a field can be non-Gaussian. In order to 
perform a thorough analysis, one therefore has to apply 
a wide range of qualitatively different tests that attack 
the problem from different directions. Different classes of 
tests, e.g., harmonic space methods such as the bispec- 
trum and trispectrum, the one-point distribution func- 
tion (e.g., skewness and kurtosis), iV-point correlation 
functions in real space, wavelet-based tests, and morpho- 
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logical measures such as those described in this paper, all 
probe different aspects of deviation from Gaussianity and 
are probably sensitive to different classes of such behav- 
ior. While tests performed in harmonic space are most 
sensitive to non-Gaussianities concentrated on character- 
istic angular scales as quantified by t, tests performed in 
real space are most sensitive to non-Gaussianities local- 
ized to specific regions on the sky. 

In this paper we present estimates of four functions, 
all related to the morphology of hot and cold areas, com- 
puted from the WMAP first-year data. Thes e four func- 
tions are the three Minkowski functionals (Minkowski 
1903) and the "skeleton" length (Novikov, Colombi, 
& Dore 2003), which all are known for their ability 
to measure non-Gaussian effects. Historically, cosmo- 
logical te sts based on topology we r e pro bably intro- 
duced by iGott. Melott. &: Dickinsonl l)1986|) in the con- 
text of the distribution of galaxies. The genus Minkowski 
functional (or equivalently, the Euler-Poincare char- 
acteristic) was first suggested as a way of detecting 
the n on-G aussianity of the CMB by iColes fc Barrowl 
l)1987|) and IColes] ((1988). Topological tests of the CMB 
anisotropics were developed further by, e.g., Gott et 
al. (1990), and finally, the full set of Minkowski func- 
tio nals was introduced into the field of CMB research 
bv [gchmalzing & Gorski ( 19981. i W ini^zki Kosowskv 
l|1998ft . and iNovikov. Feldman. fc Shandarinl l|1999j) . 
These ha ve recently been join ed by the skeleton length 
measure (jNovikov et al.ll2003D . which naturally belongs 
to the same group of statistics. One major advantage of 
these tests is that they are relatively undemanding on 
computer resources, and therefore it is feasible to apply 
these tests to the megapixel CMB data sets from modern 
experiments, such as WMAP. 

While the skeleton length has not been measured for 
the WMAP first-year data set to this date, at least three 
other groups have computed the genus Minkowski func - 
tional for this data s et, na mely, iKomatsu eTaD ll200^ . 
ICollev fc Oottl ll200l . and IParll2004D . While the for- 
mer group does not describe their algorithms in great 
depth, the latter two use pixel-by-pixel methods for com- 
puting the genus. We choose a third method based on 
the first and second derivatives of the map for finding all 
stationary points in the map. Although each of the four 
functions has an explicit analytic expectation value for 
a Gaussian field, we find it generally more convenient 
to calibrate our results with Monte Carlo simulations 
in order to assess the importance of the effects of re- 
alistic beam profiles, non-uniform noise, and partial sky 
coverage. As a by-product of these functional measure- 
ments, we also determine a number of power spectrum- 
dependent quantities, such as the number of extrema and 
the total length of the skeleton. Each of these statistics 
may be regarded as an independent test of the underlying 
power spectrum. 

It should be note d that another recent analysis by 
lHansen et alJ l)2004|) computes the local curvature dis- 
tributions of the WMAP data, using a technique very 
similar to ours. Although they do not obtain quite as 
strong detections with this curvature measure as deter- 
mined here with the genus statistic, the two sets of results 
are in excellent agreement. 

The rest of the paper is organized as follows. In |J3]wc 
review the definitions of the Minkowski functionals and 



the skeleton, while describes how to compute each 
of these functions from a pixelized map. The necessary 
preparations and the statistical methods to be applied 
to the actual WMAP data are summarized in §21 and 
[3 Finally, the results are shown in ^ before comparing 
with other results and making some final remarks in jJJ| 

2. THE ESTIMATORS 

In our analysis, we study the normalized anisotropy 
field, 

v{<P,0) = , (!) 

CO 

where the standard deviation ctq of the anisotropy field 
over the observed region 8VL of area A \^ s is defined by 

This normalization is necessary in order to minimize 
the impact of the realization-dependent power spectrum 
- we are more interested in departures from Gaussianity, 
than in possible power spectrum deviations. 

For a study of the geometric properties of the CMB 
field to be meaningful, we have to consider the shape of 
connected sets. This is done here, as usual, by studying 
excursion sets of the two-dimensional CMB-field on the 
celestial sphere, i.e., the set consisting of the parts of 
the sky having a temperature above a given threshold 
v. The main interest is in seeing how the morphological 
descriptors vary as functions of v. 

2.1. The Minkowski functionals 
A well known theorem of integral geometry l)Hadwigerl 



Il957h says that under relatively unrestrictive conditions 
(that the descriptors are rotationally and translationally 
invariant and that they preserve additivity and convex 
continuity), the morphology of a convex body in N- 
dimensional space can b e completely desc ribed by N + 1 
Minkowski functionals (Minkowski 1903). Therefore, a 
morphological description of the CMB anisotropy field 
on the two-dimensional celestial sphere can be obtained 
through three Minkowski functionals, and all possible 
morphological descriptors obeying the above conditions 
can be expressed as a linear combination of those. Thor- 
ough reviews of the two-dimensional Mink owski function- 
als on a sp here are given by, e.g., I^S chmalzing fc Gorskil 
(1998) and iWinitzki fc Kosowskvl l|1998D . 

The first two-dimensional Minkowski functional we use 
is the normalized area A of the excursion set 1Z 7 i.e., the 
area of 1Z as a fraction of the full area of the survey. Using 
the HEALPix 5 pixelization (Gorski, Hivon, & Wandelt 
1999), in which all pixels have exactly equal area, the 
area functional is easily computed by counting the num- 
ber of pixels with a field value above the threshold under 
consideration, 



dA 



Artat 
pix 



(2) 



The second two-dimensional Minkowski functional is 
the length of the border of the excursion set, 



C{v) = / di. 



(3) 



5 Available from |http: / /www.eso.org /science/healpix 
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The details of computing this measure on a pixelized map 
are provided in detail below. We follow the convention 
of normalizing the function to have unit integral value, 



J C{v)dv = 1. 



(4) 



The third Minkowski functional is the genus Q, which 
is defined as the (normalized) number of connected hot 
spots (regions above the threshold, i.e., connected parts 
of TV) minus the number of connected cold spots (re- 
gions below the threshold, i.e., connected parts of the 
whole survey area minus 1Z). However, here we approx- 
imate Q with the number of maxima plus the number 
of minima minus the number of saddle points in the ex- 
cursion set (Novikov, Schmalzing, & Mukhanov 2000), 
divided by the sum of all stationary points in the whole 
un-thrcsholded field, 



Q{y) = 



iU(^) + AUn(£) - N sad {v) 

Afmax(-Oo) + 2V min (-Oo) + iVsad(-Oo) 

2.2. The skeleton length 

The skeleton length was recently introduced a s a di- 
agnostic for Gaussianity bv lNovikov et all l|2003j) . For a 
random field on R 2 , it is defined as the zero-contour line 
of the map 



S T X Ty{T XX Ty 



T X y(T*-T*), 



(6) 



where T x and T y denote the first derivatives of the ran- 
dom field in two orthogonal directions and the -values 
denote the second derivatives. 

When applying the above expression to a field on the 
sphere, it is crucial to notice that S is a rotationally 
invariant quantity. Therefore, we may construct a lo- 
cal coordinate system at each point and compute the 
derivatives in that local system before adding them to- 
gether. In particular, we may choose our coordinate sys- 
tem to be the standard latitude-longitude system on the 
sphere, and the derivatives may be chosen to be the co- 
variant derivatives. The skeleton map S on the sphere 
may therefore be expressed by 



S — TrfT.Q^T-^ — T-gg) + T.,pg(T 2 e — T 2 



(7) 



where the semicolons as usual denote covariant deriva- 
tives. We later discuss how to compute them. 

Although Equation J7J may be difficult to interpret, 
its geometric interpretation is intuitive: the zero-contour 
lines of S are the set of lines that extend from extremum 
to extremum along lines of maximum or minimum gra- 
dient, and the set of all such lines is collectively called 
the skeleton of the field. 

The statistic we are interested in is the length of the 
skeleton of that part of the field that lies above some 
threshold v, normalized by the length of the skeleton 
of the whole un-thresholded field. Obviously this can be 
computed by the same algorithm as is used for computing 
the length Minkowski functional C, since the skeleton is 
the zero-contour of the i S-map. 

For a Gaussian field ijNovikov et al.ll2003|h the differ- 
ential skeleton length depends only on one spectral pa- 
rameter, 7, defined as 



7 



where (Tg is the variance of the map, o\ is the variance 
of the first-order derivatives, and o\ is the variance of 
the second-order derivatives. We later use this quantity 
as an independent statistic when studying the WMAP 
data. 

Actually, it has been shown that for a Gaussian field, 
the differential skeleton length is 



1 



/2tt 



[l + 0.17 7 2 (l-z> 2 )]. (9) 



3. ALGORITHMS 
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As seen from the definitions, in addition to just count- 
ing pixels for the area Minkowski functional (Equation 
0), we need two distinct algorithms for computing the 
four functions (i.e., the three Minkowski functionals and 
the length of the skeleton). For the skeleton length and 
for the length Minkowski functional it is necessary to es- 
timate the length of a contour line, and for the genus 
functional the positions of all stationary points must be 
located and the field values at these points estimated. 

Both these operations, as well as the determination of 
the skeleton, depend on derivatives of the field. Such 
derivatives are well-behaved only if the map itself is 
smooth (that is, not dominated by pixel noise), and we 
therefore consider it prudent to filter each map before 
applying our algorithms. In addition, simple experi- 
ments show that robust estimation of the contour line 
lengths requires a high pixel resolution compared to the 
smoothing scale: adequate smoothing is then an essen- 
tial step in the analysis. A brief discussion of this and 
related topics is given in ® As pointed out by, e.g., 
iWinitzki fc Kosowskvl l)1998j) . numerical techniques for 
estimating the length of a contour line or positions of 
stationary points on real CMB maps with noise, pix- 
elization, masks, etc., have to be chosen with care. Our 
experiments show that the algorithms employed here are 
sufficiently insensitive to such effects and give accurate 
results. 

3.1. Measuring the length of an iso-contour line 

Our starting point is thus a smoothed, pixelized 
map in which we know the field values at the pixel 
centers. We want to trace the underlying contour 
lines between those pixel centers and estimate their 
lengths. The method for doing so is based on lin- 
ear interpolation and is almost identical t o the meth- 
ods adopted by. e.g.. jShandarin et alJ lj2002j) and 
iNovikov. Feldman. fc Shandarinl l)1999|) . However, we 
also use this method for locating the stationary points; 
therefore, we review the algorithm in detail here. 

We start by constructing a secondary set of pixels, de- 
fined by letting the original pixel centers be at the corners 
of the new, secondary pixels. The secondary pixels are 
then checked to determine those that are crossed by a 
contour line. These are obviously the secondary pixels 
where at least one vertex has a lower value and at least 
one vertex has a higher value than the contour line. This 
process is illustrated in the left-hand panel of Figure ^ 

Once such a "crossed" secondary pixel has been found, 
it is easy to determine through which edges the contour 
line passes; a contour line crosses an edge if one vertex 
has a field value larger than the contour value and the 
other vertex has a smaller value. The point on the edge 
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a) b) 

Fig. 1. — Tracing contour lines and localizing stationary points by means of linear interpolation, (a) Temperature values Ty are known 
only at the centers of the HEALPix pixels, and linear interpolation is therefore used to approximate the true contour line. The true contour 
line is shown as a dashed curve in this figure, and the linear approximation as solid line segments. In this process it is useful to construct 
a set of secondary pixels defined by letting the centers of the HEALPix pixels define their vertices. In order to trace the contour line, one 
then has to check all such secondary pixels, searching for edges for which the two vertices have different signs relative to the contour line. 
Once such an edg e has been located, the point of intersection (marked as p; in the figure) is approximated by linear interpolation, according 
to Equation I1UI . (b) Localizing stationary points is similar to tracing the contour lines. In this case one focuses on the zero-contours of 
the first derivative maps and searches for secondary pixels in which both der ivat ives are zero. Once such a pixel is found, the position of 
intersection is approximated by linear interpolation, i.e., by solving Equation 1171 for s and t. If < s,t < 1, then the point of intersection 
lies within the secondary pixel under consideration. 



where the contour line crosses is estimated by a linear 
approximation, 



AT, - AT C , 



AT 2 - ATi 



ni 



ATi - AT cl 



AT 2 - ATi 



n 2 . (10) 



Here n cont (given as an unnormalized vector from the 
center of the sphere) is the point at which the contour 
line crosses the secondary pixel edge, ATi and AT 2 are 
the field values at the two vertices with corresponding 
positions ni and n 2 , while AT cont is the field value on 
the contour line. 

After locating both points on the edges of the sec- 
ondary pixel where the contour line enters and exits the 
pixel, the angular length of the line segment within that 
pixel is estimated by computing the dot product of the 
two vectors, 



5£ — arccos(nJ Q 



n 2 ) 



(ii) 



where n* ont are the normalized unit vectors found in 
Equation jlOJ. 

In those cases in which more than two edges of a pixel 
are crossed by a contour line, we assume that two contour 
lines cross each other inside the pixel. The justification 
for this assumption is simply that we are interested only 
in the total length of the contour line, not its shape, and 
the possible error introduced in the total length by this 
assumption is obviously very small because of the small 
pixel size. 

3.2. Locating the stationary points of a random field on 

a sphere 

In order to compute the genus Minkowski functional, 
we need to locate all stationary points - both minima, 
maxima, and saddle points. While extrema alone may be 
found by searching for pixels that have neighbors with all 
higher or lower temperature values (as is implemented in 



the HEALPix hotspot facility), this method cannot be 
used for locating saddle points. Therefore, a different ap- 
proach is adopted in this work, namely, a method based 
on the covariant derivatives of the random field. 

A stationary point has by definition vanishing first par- 
tial derivatives computed in two orthogonal directions. 
The particular type of stationary point is determined by 
the eigenvalues of the Hessian matrix: if all eigenvalues 
are positive, the point is a minimum; if all are negative, 
the point is a maximum; and if there are both positive 
and negative eigenvalues, it is a saddle point. To find the 
positions and types of all stationary points, we therefore 
need to be able to efficiently compute the first and second 
derivatives of the temperature field. We use the covariant 
derivatives in the standard polar coordinate system, 



AT fl 



1 9 AT 

sine? def) 
dAT 
89 ' 



(12) 
(13) 



In this paper, these derivatives are only used to deter- 
mine the skeleto n and to locate the s t ationa ry points, 
but, as shown bv lSchmalzing fc Gorskil (^998), it is pos- 
sible to express all the Minkowski functionals in terms of 
the field itself and its first- and second-order derivatives. 

It is most convenient not to compute the covariant 
derivatives of the temperature anisotropy field directly 
from the pixel values, but rather to do it from the ex- 
pansion of the anisotropy field in spherical harmonics, 

/max I 

AT{6, </>) = g 12 ai ™ Yl ™( 9 > & ( 14 ) 

1=0 m=-l 

'max I 

= 12 12 a lm X lm P lm {cos9)e lm \ (15) 



2=0 ro=-i 
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Galactic longitude, I Galactic longitude, / 



a) b) 

Fig. 2.— A 10° X 10° patch from the CMB sky as measured by WMAP, with (a) dT/d6 = marked in blue and dT/d(f> = marked in 
red, and (f>) the skeleton of the field marked in yellow. In both figures red filled circles indicate maxima, blue filled circles indicate minima 
and green crosses indicate saddle points. 



where the partial derivatives are found term by term. 

Here X lm = {[{21 + 1)/4tt] [(/ - m)!/(Z + m)!]} 1/2 is a 
normalization constant. 

By using the identity for associated Legendre polyno- 
mials, 



dPT 

d.r 



1 



1 — x 



{l + l)(! + m) 
21 + 1 



P, 



1(1 -m + l) 



i-i 



Pi' 



21 + 1 ' 1+1 
(16) 

where x = cos 9, one finds that the covariant derivatives 
are given by 



AT : „ 



AT :8 = 



E E (ma lm )Y lr 

1 = rn=-l 



sin 6 



E 



l—i 

E 

n=-l+l 



-1 



I 



E 

1=0 



E 



(i - - m) 
21 - 1 ' 

(i + 2)(i + m + 1). 
21 + 3 



A"/ 1 



A, 



-Ol + l ! 



A,, 



These derivatives may then be computed very efficiently 
using the HEALPix routines map2alm and alm2map. 

In principle, one could go through this process once 
more to find equivalent expressions for the second deriva- 
tives. However, in practice one obtains numerically more 
stable results by making new pixelized maps of the first 
derivatives of the original map, expanding those in spher- 
ical harmonics using HEALPix routines, and then ap- 
plying the above formulae to them, rather than using 
second-order expressions on the temperature map itself. 
The cost for this extra stability comes in the form of 
two extra spherical harmonics transforms, but the extra 
CPU time required is acceptable. By repeated appli- 
cation of the above formulae we thus find all first and 
second derivatives of the anisotropy field in all pixel cen- 
ters. 

The next problem is to find the points of intersection 
of the lines AT ; ^ = and AT-e = 0, since the stationary 
points must be at these intersections. This is done with 



an algorithm similar to that described for tracing the iso- 
contour lines; once again we focus on the set of secondary 
pixels, but this time we search for secondary pixels being 
crossed by both "iso-contour lines" AT.^ = and AT.g = 
0. 

Once such a secondary pixel is found where the two 
zero-derivative curves cross its edges, it must be deter- 
mined whether the two lines actually cross each other 
inside that pixel. We find the position of all four edge 
crossings (two for each contour line) as for Equation (|1UII , 
project these positions into a local, two-dimensional coor- 
dinate system centered on the pixel, and solve the system 



spi+(i-s)pt t =tpi+(i-t)pt 



of linear equations. Here the four p* 



in/out 



(17) 
are four 

two-dimensional vectors representing each of the four 
crossing-points of the edge of the pixel. 

If < s, t < 1, the two lines do intersect each other 
inside the pixel, and the point of intersection is estimated 
by substituting either s or t into the respective side of 
the above equation. Figure ^b) illustrates this method. 

After locating such a stationary point, one must then 
estimate both the field value (for thresholding) and the 
second derivatives at that point. For this purpose we 
adopt a weighting method that takes into account the 
values at all four vertices. First, the angular distance 
between the point of intersection and the four vertices is 
computed, 

m — arccos(p • pi). (18) 

The relative weights, u)j, are then evaluated as the small- 
est angular distance divided by Ui, 



(19) 



Finally, the true weights, lUj, are found by normalizing 
the sum to unity, ^]w t = 1. 
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Fig. 3. — Masks used in the computations of the Minkowski 
functionals and the skeleton for FWHMs between 0?53 and 4? 26. 
The masked regions at high latitudes are extended "semi-point" 
sources excluded by the KpO mask. 



The field value and the second derivatives at the point 
of intersection may now be estimated as weighted sums 
over the vertices, 

4 

AT = Y^ Wi AT 1 (20) 

4 

AT ;jh »< Ar - ( 21 ) 

The last step in this process is to determine the type 
of stationary point. This is done from the characteristic 
equation, 

\ 2 -(AT.0g+AT.^)\+{AT.g 9 AT.^-AT^) = 0. (22) 

The stationary point is a maximum if the solutions of 
the characteristic equation Ai,A2 < 0, a minimum if 
Ai,A2 > and a saddle point if Ai and A2 have differ- 
ent signs (assuming that the point is non-degenerate, i.e., 
that \i which is always true for real- world fields). 

The procedure is illustrated in Figure J2J. A contour 
map of a patch of the CMB sky is shown in the left 
panel. The APg = contours are the blue lines and 
the AT ; = contours are the red lines. The stationary 
points (maxima: red filled circles; minima: blue filled 
circles] saddle points: green crosses) are seen to be at 
the intersections of the blue and red lines. In the right- 
hand panel the skeleton of the field is overlaid on the 
same map as a yellow contour line. 

4. THE WMAP DATA AND THE SIMULATIONS 

The aim of this paper is to measure the three 
Minkowski functionals and the length of the skeleton as 
a function of the threshold height for the WMAP first- 
year maps and to compare these to the same quanti- 
ties estimated from a Monte Carlo set of 5000 simulated 
Gaussian sky maps. The Monte Carlo ensemble is based 
on the best-fit WMAP power spectrum with a running 
index, which may be down-loaded from the LAMBDA 6 
Web site. 

6 http://lambda.gsfc.nasa.gov 



Our study is based on the co-added WMAP map (Hin- 
shaw et al. 2003), and in order to produce the best possi- 
ble approximations to this map, the simulations are pro- 
cessed in exactly the same way as the WMAP map. In 
short, each realization is constructed as a weighted sum 
over each of the cosmologically important WMAP chan- 
nels (Ql, Q2, VI, V2, Wl, W2, W3, and W4), taking 
into account t he particular beam and noise properties of 
each channel {Spcrgci et al. 2003J). We use foreground- 
c orrected maps produced by the template fitting method 
of lBennett et al.l l)2003b|) to suppress foreground contam- 
ination. 

The publicly-released maps are available at the 
LAMBDA Web-site. However, a problem was recently 
discovered with these maps (see the LAMBDA web 
page 7 , where a more thorough description is given); the 
templates used in this process were convolved with a too 
low / max (the maximum multipole component) , resulting 
in ringing around strong point sources. For this reason 
we have also generated our own templates and performed 
the analysis for these maps as well. The differences be- 
tween the two sets of results are insignificant on the scales 
we consider; therefore, we choose to report only the re- 
sults for the official maps in what follows. 

For both the WMAP data and the simulated maps, a 
number of data processing steps are performed to min- 
imize the impact of residual foregrounds, the applied 
Galaxy cut, pixel noise, etc. First, the WMAP KpO mask 
with point sources included (i.e., we do not exclude the 
700 point sources present in the mask specified by the 
WMAP team) is defined as the base mask. The decision 
not to exclude the 700 point sources is based on two fac- 
tors. First, we are mainly interested in scales larger than 
2°, and for such large scales the impact of point sources 
is small (see, e.g., Spergel et al. 2003). Later in the paper 
the impact of point sources is explicitly tested by apply- 
ing a median filter to the map. Second, since the maps 
are smoothed, we also have to expand the mask accord- 
ingly. If the 700 point sources actually were excluded by 
the base mask, the expanded mask would accept only a 
very small fraction of the sky. It is therefore better to 
include the point sources in the analysis, and then later 
on determine from the results themselves whether this is 
justified or not. To test for any effects of the exact shape 
of the mask, we also perform our analysis using a base 
mask with the additional constraint \b\ < 30°. 

The next processing step is to remove the best-fit 
mono- and dipoles from the maps, with coefficients com- 
puted from the accepted region only. Then, the masked 
region is nullified, and the mask boundary is apodized, 
before the maps are smoothed with a Gaussian beam of 
given FWHM (which will be varied, see below). This 
smoothing operation consists of a spherical harmonics 
transform of the original map, followed by a multiplica- 
tion with the Legendre transform of the Gaussian beam. 
Then the smoothed map is constructed through an in- 
verse spherical harmonics transform. At this point we 
keep the final multipole expansion components, a; m , for 
computing the derivatives (see © at a later stage. 

The next step is to expand the excluded regions of 
the mask in all directions as a safeguard against bor- 
der effects. This is done according to the procedures of 

7 http:/ /lambda. gsfc.nasa.gov/product/map/IMaps_cleaned.cfm 
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iSchmalzing &: Gorskil l)1998|k the mask map consisting 
of zeros and ones is convolved with a Gaussian beam 
of the required FWHM, and the new, expanded mask 
is then determined by including all pixels with values 
higher than some given threshold, for instance, 0.95 or 
0.99. For maximum safety we perform this operation 
twice, each time with a limiting threshold of 0.99. The 
resulting masks are shown in Figure |3 Finally, the map 
is normalized according to Equation ^ including pixels 
in the accepted region only. 

The analysis is carried out for a large number of 
smoothing FWHMs in order to probe different angu- 
lar scales. The beam widths selected were 0?53, 0?64, 
0?85, 1?28, 1?70, 2?13, 2?55, 2?98, 3?40, 3?83, and 4?26 
FWHM 8 . Note that this smoothing is applied directly to 
the observed maps; therefore, the Gaussian smoothing 
beam is in addition to the experimental WMAP beams. 
In total, the narrowest of the beams (0?53 FWHM) is 
sensitive to multipoles up to about I w 600; therefore, 
both noise and point sources are expected to play a ma- 
jor role at this and s maller scales. Here it is worth notic- 
ing that lParkl l)2004|) considered even smaller scales, since 
the maps were not explicitly smoothed at all, except for 
a small effect introduced by a stereo-graphical projec- 
tion. However, in the subsequent analysis, the 700 point 
sources resolved by WMAP were masked out, and it was 
therefore possible to produce results on smaller scales 
than presented in this paper. 

The maximum multipole component, £ max , in each run 
is chosen to match the corresponding FWHM. In partic- 
ular, the values of l max we choose are 650, 600, 550, 
and 450 for FWHMs between 0?53 and 1?28 and 350 
for the larger beams. The HEALPix resolution parame- 
ter used internally in the computations is N s ^ = 1024. 
By increasing A^ide from 512 to 1024, the total skeleton 
length is increased by a few percent for the narrowest 
beam. However, by increasing Aside from 1024 to 2048, 
the total skeleton length is increased only by a negligi- 
ble amount. Thus, A s id c = 1024 is the lowest resolution 
able to support all our analyses, although lower resolu- 
tions could have been chosen individually for the wider 
beams. 

The four statistics are estimated for 200 values of the 
threshold height between — 4ero and 4cto, but only a sub- 
set of these values are used in each of the subsequent 
analyses. Each statistic is evaluated independently on 
the northern and southern Galactic hemispheres. Inter- 
esting hemisphere effects have been reported by several 
authors (Eriksen et al. 2004; Park 2004), and we focus 
our tests accordingly. 



5. QUANTIFYING THE DEGREE OF AGREEMENT 
BETWEEN SIMULATIONS AND OBSERVATIONS 

In order to quantify the degree of agreement between 
the simulations and the observations, we adopt a diago- 

8 After having completed the computations, we discovered an 
error in a computer code that effectively multiplied all smoothing 
FWHMs by V8 log 2 = 2.35. Thus, the FWHMs that originally 
had been chosen to be multiples of 1° were in reality multiples of 
0?43, causing the somewhat unnatural looking smoothing scales in 
this paper. 



nal x 2 statistic of the form 



X 



E 

£>£(-3,3) 



m - (m) 



n 2 



a{v) 



(23) 



where v is the threshold level, (/(£')) is the average of 
the simulations, and cr(v) is the standard deviation of the 
measure under consideration. This statistic is computed 
for both the simulations and the observed functions, and 
the fraction of simulations with x 2 value lower than the 
X 2 determined from the real WMAP map is connected 
to the confidence level at which to either accept or reject 
the null-hypothesis - that the observed field is a natural 
member of the Gaussian Monte Carlo ensemble. 

Note that we choose a diagonal x 2 statistic for the func- 
tions in this paper. In our experience, the ordinary % 2 
statistic, which takes into account bin-to-bin correlations 
through the covariance matrix, behaves rather badly if 
the functions under study are binned with a narrow bin 
size, as is the case here. If neighboring bins are strongly 
correlated, the covariance matrix converges very slowly, 
and in the study here, 5000 simulations would not be 
sufficient to produce robust results. Furthermore, when 
the bin-to-bin correlations are strong, the off-diagonal 
terms become more and more dominant, and the shape 
of the inverse covariance matrix resembles that of a Mex- 
ican hat. A Mexican hat matrix is the same as a high- 
pass filter. Therefore, if the bin-to-bin correlations are 
strong, the covariance matrix % 2 statistic is more sensi- 
tive to fluctuations in the functionals, rather than abso- 
lute deviations. In some cases, one actually finds that 
the covariance matrix x 2 statistic accepts functions that 
should obviously be rejected by eye. 

However, if the proposed model for the simulations is 
an adequate description of the real data, then no statistic 
should be able to distinguish between the simulations and 
the observations. The choice of statistic is therefore only 
a question of what one wants to measure. 

For the genus we adopt the machinery of Park ( 2004) 
and perform a parametric fit of the form 

Getty) = A(y - AD)e-^- A ^ 2 , (24) 
where A and Av are free parameters. For a Gaussian 
field it is expected that AD — 0, while A is a nor- 
malization constant depending on the power spectrum. 
This fit is performed over three different ranges in D, 
namely, D = (-3,3), (-2.5,-0.2), and (0.2,2.5). The 
corresponding amplitudes arc named A, A_, and A + 
and indicate, respectively, the best-fit total amplitude 
and negative and positive amplitudes. We also define an 
asymmetry parameter by 

Ag= A -~ A A + . (25) 
A- + A + y ' 

For a Gaussian field there should be no asymmetry in 
the shape of the genus, and consequently Ag should also 
be zero. No te that this definition of Ag differs slightly 
from that of lParkl l|2004[) . 

The fit of Get{v) to the data is implemented as a non- 
linear search (using NAG 9 routines), which minimizes 
the functional 



Go\>s{v) - Gs.t{v) 



(26) 



9 Numerical Algorithms Group, http://www.nag.co. uk 
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Fig. 4. — Minkowski functionals and differential skeleton measured from the KpO sky, smoothed with a 1?28 FWHM Gaussian beam. 
Gray bands indicate 1, 2, and 3<r bands, and the solid line indicates the ensemble average. Filled circles show the observed functions. 



where w(y) is a weight factor. In this paper we choose in- 
verse variance weighting, w{v) — , where Og{v) 
is the simulated variance of the particular genus bin. 
However, almost all results are practically independent 
of this particular choice, with the notable exception of 
the total range amplitude A for realizations with a large 
asymmetry parameter Ag. In those cases, one may find 
that the amplitude depends, but only weakly, on the 
weighting scheme. Moreover, we are mostly concerned 
with the negative and positive threshold amplitudes and 
the asymmetry parameter itself, and these quantities are 
very robust against different weighting schemes. 

6. RESULTS 

The algorithms described in J|l result in a number of 
by-products, and several of these may serve as useful 
consistency checks. In particular, the validity of the as- 
sumed power spectrum may be tested by e.g., counting 
the number of stationary points, or by measuring the 
total length of the un-thresholded skeleton. This is the 
topic of the next subsection, while the results from the 
actual non-Gaussianity analysis are presented in the fol- 
lowing subsection. 

In the following, the results from the measurements 
are presented as a function of the smoothing scale. One 
should therefore remember that each quantity is a con- 
tinuous function of that smoothing scale, and the results 



at two different scales cannot be considered independent. 

6.1. Power spectrum consistency statistics 

In Tabled the spectral index 7, the skeleton length, 
and the stationary point count results for the un- 
thresholded map using KpO as the base mask are shown, 
for both the full sky and for the northern and south- 
ern hemispheres separately. The numbers in parentheses 
show the percentage of simulations with a lower value 
than that found in the WMAP sky. The value 100% is 
reserved for the special case in which all 5000 simulations 
have a lower value. Note also that the percentages refer 
directly to the cumulative distribution functions and not 
to the significance level as such; both 2.5% and 97.5% 
indicate a 2 a effect. 

On the smallest scales (0?53 and 0?64 FWHMs) the 
spectral parameter 7 is high, but well within the accept- 
able range. This is also the case for the total length of the 
skeleton. However, the numbers of stationary points are 
clearly not acceptable - only two simulations out of 5000 
realizations have as many saddle points as the WMAP 
map at 0?53 FWHM, and the overall stationary point 
counts lie at around 3 a compared with the distribution 
of the 5000 simulated results. However, this is not a 
very unsettling result, considering that our choice of base 
mask does not discard the 700 known point sources. It 
should be expected that point sources cause deviations 
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TABLE 1 

Power spectrum dependent measurements, KpO mask 



FWHM 



Length 



N min 



Full sky 



0?53 
0°64 
0?85 
1?28 
1?70 
2? 13 
2? 55 
2? 98 
3? 40 
3? 83 
4? 26 



0.497 (95.2%) 
0.505 (93.9%) 
0.503 (93.6%) 
0.478 (93.4%) 
0.468 (93.4%) 
0.465 (93.9%) 
0.470 (95.0%) 
0.478 (96.4%) 
0.488 (97.2%) 
0.496 (97.7%) 
0.505 (98.0%) 



1870. 
1584. 
1267. 
904.6 
677.8 
530.1 
429.7 
354.2 
297.0 
254.4 
221.2 



(84.3%) 
(93.4%) 
(93.0%) 
(67.9%) 
(58.6%) 
(84.0%) 
(97.6%) 
(85.7%) 
(56.0%) 
(52.0%) 
(54.1%) 



11258 
8558 
5723 
2939 
1648 
1042 
691 
494 
366 
285 
216 



(98.4%) 
(99.5%) 
(91.8%) 
(25.9%) 
( 6.5%) 
(49.8%) 
(48.2%) 
(73.4%) 
(82.5%) 
(94.9%) 
(80.5%) 



11330 
8515 
5701 
3019 
1716 
1063 
705 
499 
348 
264 
204 



(99.9%) 
(96.2%) 
(81.7%) 
(97.5%) 
(88.9%) 
(85.9%) 
(80.7%) 
(84.7%) 
(28.5%) 
(35.3%) 
(34.0%) 



22 565 
17 049 
11422 
5958 
3352 
2100 
1395 
986 
719 
554 
422 



(99.9%) 
(99.3%) 
(92.2%) 
(80.5%) 
(32.9%) 
(71.3%) 
(70.0%) 
(78.7%) 
(72.0%) 
(87.7%) 
(70.1%) 



Northern hemisphere 



0?85 
L°28 
1?70 
2? 13 
2? 55 
2? 98 
3? 40 
3? 83 
4? 26 



0.535 (98.7%) 
0.517 (98.9%) 
0.511 (98.9%) 
0.514 (99.0%) 
0.526 (99.3%) 
0.537 (99.5%) 
0.548 (99.6%) 
0.555 (99.3%) 
0.559 (98.9%) 



634.5 
454.8 
334.0 
268.2 
217.2 
177.3 
149.1 
128.6 
113.0 



(70.3%) 
(60.6%) 
(23.9%) 
(83.6%) 
(89.3%) 
(35.5%) 
(22.6%) 
(33.3%) 
(61.7%) 



2851 (62.4%) 
1493 (57.6%) 
811 ( 1.3%) 
524 (41.9%) 
339 (16.3%) 
245 (45.0%) 
181 (52.7%) 
143 (79.5%) 
112 (80.6%) 



2854 (65.9%) 
1537 (98.6%) 
875 (92.7%) 
545 (90.5%) 
360 (80.9%) 
243 (37.0%) 
170 (10.2%) 
138 (57.5%) 
105 (42.6%) 



5694 (60.4%) 
3023 (91.2%) 
1681 (27.2%) 
1068 (75.9%) 
702 (56.9%) 
490 (47.6%) 
356 (40.0%) 
284 (83.7%) 
218 (72.6%) 



Southern hemisphere 



0?85 

r?28 

1?70 
2? 13 
2? 55 
2? 98 
3? 40 
3? 83 
4? 26 



0.492 (27.9%) 
0.463 (25.5%) 
0.449 (25.0%) 
0.444 (25.0%) 
0.445 (26.9%) 
0.452 (31.2%) 
0.460 (36.4%) 
0.471 (42.9%) 
0.485 (50.7%) 



632.5 
449.8 
337.7 
261.9 
212.6 
177.2 
147.6 
125.7 
108.1 



(94.3%) 
(64.4%) 
(83.2%) 
(65.5%) 
(91.9%) 
(97.4%) 
(77.5%) 
(63.1%) 
(41.1%) 



2871 (94.2%) 
1453 (19.8%) 
834 (45.0%) 
517 (52.4%) 
346 (64.1%) 
252 (89.3%) 
187 (91.6%) 
139 (82.3%) 
104 (57.9%) 



2851 (83.0%) 
1477 (59.9%) 
832 (38.8%) 
517 (51.6%) 
345 (58.6%) 
259 (97.9%) 
175 (48.3%) 
127 (23.7%) 
98 (23.9%) 



5719 (94.0%) 
2937 (44.7%) 
1661 (34.2%) 
1029 (45.8%) 
688 (59.6%) 
501 (92.3%) 
365 (86.5%) 
268 (65.5%) 
204 (50.2%) 



Note. — The table gives the spectral parameter 7, the length of the skeleton (in 
radians), the number of minima, maxima, and saddle points as a function of beam 
FWHM for the un-thresholded co-added WMAP map using the KpO mask. The numbers 
in parentheses are the percentage of the 5000 simulations with a lower value than that 
found in the WMAP map. The results are for the full map, and for the northern and 
southern hemispheres separately. 



TABLE 2 

Power spectrum dependent measurements, KpO mask using only |6| > 30° 



FWHM 



Length 



Northern hemisphere 



2? 13 
2? 55 
2? 98 
3? 40 
3? 83 
4? 26 



0.538 (98.7%) 
0.549 (99.1%) 
0.561 (99.1%) 
0.565 (98.1%) 
0.566 (95.8%) 
0.569 (92.9%) 



167.7 (57.5%) 
134.9 (75.4%) 
109.6 (38.7%) 
92.2 (54.3%) 
78.9 (69.6%) 
69.4 (92.3%) 



327 (35.7%) 
215 (34.3%) 
152 (47.0%) 
117 (82.7%) 
89 (84.1%) 
74 (97.3%) 



339 (75.5%) 
219 (52.5%) 
153 (53.3%) 
106 (23.8%) 
86 (70.4%) 
67 (71.4%) 



664 (55.2%) 
437 (52.0%) 
311 (71.7%) 
223 (60.0%) 
173 (80.5%) 
140 (94.3%) 



Southern hemisphere 



2? 13 
2? 55 
2? 98 
3? 40 
3? 83 
4? 26 



0.471 (38.5%) 
0.476 (40.8%) 
0.480 (40.2%) 
0.488 (42.7%) 
0.498 (45.0%) 
0.505 (45.4%) 



166.1 (52.7%) 
134.9 (83.2%) 
113.6 (98.2%) 
95.0 (91.2%) 
80.7 (78.7%) 
68.9 (44.5%) 



328 (48.9%) 
227 (84.7%) 
168 (97.6%) 
126 (98.2%) 
87 (67.9%) 
70 (79.8%) 



328 (48.0%) 
225 (78.5%) 
172 (99.3%) 
117 (78.3%) 
85 (53.5%) 
63 (31.5%) 



651 (38.3%) 
453 (89.8%) 
336 (99.3%) 
244 (97.2%) 
180 (89.6%) 
137 (78.4%) 



Note. — Same as Table but where all areas with |6| < 30° have been 
included in the mask, i.e., excluded from the analysis. 
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on these scales. In fact, these measurements are only in- 
tended to set a meaningful lower limit on the smoothing 
FWHM to use in the following analyses. 

At the other extreme, there are only a few hundred 
stationary points in the full sky at the 4? 26 scale. Since 
the method for computing the genus is based on count- 
ing stationary points, we must expect artifacts at larger 
scales; therefore, we choose not to extend our analysis to 
this regime. 

On scales larger than 1°, all full-sky numbers are in 
fairly good agreement with the model, although the spec- 
tral parameter 7 is somewhat large (> 93% of the simu- 
lations). However, the situation is much more interesting 
when looking at the northern hemisphere separately. In 
this case the WMAP spectral parameter 7 is higher than 
98.9% of the simulation on all scales, reaching the value 
of 0.548 at 3?40, which is high at the 99.6% level. Ex- 
cept for the spectral parameter, all measurements are 
perfectly acceptable. No sign of discrepancy is found in 
the southern hemisphere. 

In Table |21 we show the same statistics for a number 
of smoothing scales in the northern and southern hemi- 
spheres when all areas with \b\ < 30° are added to the 
base mask. Extending the excluded region only alters the 
results by small amounts in the northern hemisphere, but 
in the southern hemisphere the stationary point counts 
are relatively high for FWHMs around 3°. 

6.2. The Minkowski junctionals and the skeleton length 

The three Minkowski functionals and the length of the 
skeleton, computed from maps smoothed with a 1?28 
Gaussian beam (KpO base mask), are shown in Figure 
0| The WMAP data are plotted with black filled circles, 
while the median computed from the ensemble of 5000 
Monte Carlo simulations is shown as a thin solid line. 
The 1, 2, and 3(7 confidence bands around the median 
are shown as gray bands. However, in all cases except 
the genus, the confidence regions are so narrow that it is 
virtually impossible to distinguish between the different 
elements. 

For this reason, we subtract in what follows the median 
from each bin before plotting the functions, and such re- 
sults are shown for the full-sky measurements in Figure 
03 In these plots there are several points worth noticing: 
First, we see that the genus confidence bands are asym- 
metric in the tails. This is because of the inability of our 
method to deal with masks; i.e., on a cut sky, the number 
of maxima plus the number of minima minus the number 
of saddle-points is likely to be non-zero, and this trans- 
lates in our method into a non-zero genus at v = —00. 
On the other hand, there are no stationary points at very 
high thresholds (when the remaining area of the sky is 
very small); therefore, the genus will always be zero here. 
However, this asymmetry is only a very minor problem 
for thresholds between — 3oo and 3ctoj an d the effect is 
in any case calibrated by simulations. 

Second, the plots show that the properties of the full- 
sky field on scales up to 3° are in fairly good agreement 
with the model. No single function shows excessive de- 
viations, but the deviations are also not too small com- 
pared to the confidence bands, as could easily be the case 
if the power spectrum was flawed. However, the genus 
curve does seem to have a somewhat large amplitude, 
and more so at large scales. At 3? 40 FWHM the genus 



TABLE 3 

Minkowski functional and length of 

SKELETON \ 2 RESULTS, KpO MASK 



FWHM Area Length Genus Skeleton 



Full sky 



0?53 


49.4% 


67.8% 


86.2% 


23.4% 


0?64 


43.9% 


65.2% 


88.3% 


28.6% 


0?85 


35.2% 


56.8% 


89.4% 


35.4% 


1?28 


26.9% 


44.3% 


89.3% 


29.3% 


1?70 


23.5% 


34.8% 


86.0% 


28.4% 


2? 13 


21.5% 


29.1% 


79.5% 


24.8% 


2? 55 


19.9% 


19.1% 


77.5% 


20.8% 


2? 98 


16.2% 


13.8% 


69.3% 


16.6% 


3? 40 


12.5% 


10.8% 


95.9% 


11.3% 


3? 83 


8.2% 


7.6% 


94.3% 


8.1% 


4? 26* 


4.3% 


3.1% 


80.5% 


7.8% 


Northern hemisphere 


0?85 


12.5% 


34.5% 


95.9% 


14.4% 


1?28 


15.4% 


32.5% 


95.3% 


15.1% 


1?70 


14.1% 


24.6% 


97.2% 


18.9% 


2? 13 


9.6% 


26.7% 


95.7% 


13.8% 


2? 55 


9.4% 


30.6% 


98.4% 


18.4% 


2? 98 


13.1% 


34.6% 


98.4% 


23.0% 


3? 40 


17.7% 


35.3% 


99.5% 


20.3% 


3? 83 


20.3% 


29.6% 


94.9% 


32.2% 


4? 26* 


28.3% 


29.4% 


77.3% 


40.1% 


Southern hemisphere 


0?85 


33.5% 


54.9% 


27.7% 


33.6% 


1?28 


53.2% 


70.9% 


37.7% 


58.4% 


1?70 


59.2% 


64.0% 


43.0% 


67.0% 


2? 13 


58.6% 


65.7% 


28.8% 


55.9% 


2? 55 


58.6% 


71.4% 


19.4% 


54.6% 


2? 98 


54.2% 


70.5% 


85.2% 


49.5% 


3? 40 


51.4% 


66.6% 


57.6% 


43.9% 


3° 83 


45.0% 


56.3% 


21.8% 


35.5% 


4? 26* 


39.6% 


40.4% 


9.6% 


32.2% 



Note. — Area functional, length func- 
tional, and genus results from diagonal x 2 
tests, as computed from values of the thresh- 
old v between — 3 00 and 3cro. The numbers 
indicate the percentage of simulated realiza- 
tions with a x 2 value lower than that for the 
co-added WMAP map using the KpO mask. 

*X 2 computed only for thresholds between 
—2.5 cro and 2.5 co- 



kes clearly outside the 2 a band. 

In order to quantify the degree of deviation in each 
case, we use the diagonal x 2 statistic described in <J5J 
The results from these measurements are shown in Table 
(for the KpO base mask) and Table 0] (for the larger 
base mask). The numbers in these tables indicate the 
percentage of simulations with a lower \ 2 value than 
that of the observed WMAP data. Once again we see 
that the full-sky numbers are in excellent agreement with 
the Gaussian theoretical model, except for the genus at 
large (3°-4°) scales, which has a x 2 that is high at the 
96% level. 

This pattern is clearly stronger in the northern hemi- 
sphere. Here the genus x 2 generally is very large, with 
a pronounced peak at 3?40. At that scale only 0.5% of 
the simulations have a larger x 2 ■ The other functionals 
are all perfectly acceptable. In the southern hemisphere, 
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Fig. 5. — Results from the Minkowski functional and skeleton measurements, with the median subtracted from each bin. Without such a 
subtraction it is not possible to distinguish between a function that follows the median and one that lies in the 2 a region, as a result of the 
extremely narrow confidence bands - usually only a few percent wide. The gray bands indicate 1 and 2 a confidence regions as computed 
from the simulations. 



even the genus is acceptable, and no signs of discrepancy 
between the observed sky and the simulations are found 
there. We see that these results are not particularly de- 
pendent on the choice of base mask, except that with the 
larger mask, the genus x 2 peaks at 2? 98 FWHM in the 
northern hemisphere, rather than at 3?40, and that it is 
curiously low at 2?55 in the southern hemisphere (at the 
1% level). However, it is difficult to attach much signifi- 



cance to the latter result, since it is highly unstable with 
respect to both Galactic cut and smoothing scale. 

The x 2 statistic only measures the overall level of devi- 
ation between a given function and the theoretical mean. 
It is possible that other statistics could be more sen- 
sitive to particular features. To study the behavior of 
the genus in more detail, we therefore fit the para metric 
function given by Equation (|24|l to each realization ijParkl 
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Fig. 6. — Large-scale genus measurements. The northern hemisphere genus functional at 3? 40 FWHM has a x 2 value that is larger than 
that found in 99.5% of the simulations, and only 1 out of 5000 simulations has a larger negative threshold amplitude. Conversely, the 
southern hemisphere genus functional at 2?98 FWHM has a very small amplitude, and its asymmetry parameter is extreme at the 99% 
level. 



TABLE 4 

Minkowski functional and skeleton x 2 

RESULTS, KpO MASK USING ONLY [6| > 30° 



FWHM Area Length Genus Skeleton 



Northern hemisphere 



2? 13 


22.5% 


48.0% 


97.6% 


22.6% 


2? 55 


11.3% 


49.8% 


99.2% 


20.2% 


2° 98 


9.2% 


49.8% 


99.9% 


33.2% 


3? 40 


10.6% 


46.0% 


94.5% 


29.5% 


3? 83 


12.3% 


40.2% 


62.2% 


24.5% 


4? 26* 


12.1% 


27.5% 


42.0% 


16.5% 


Southern hemisphere 


2? 13 


31.5% 


18.7% 


16.6% 


35.9% 


2? 55 


36.2% 


12.7% 


1.0% 


36.8% 


2? 98 


31.5% 


7.2% 


34.1% 


36.3% 


3? 40 


34.6% 


13.6% 


17.1% 


25.0% 


3? 83 


33.5% 


16.2% 


96.8% 


19.2% 


4? 26* 


32.0% 


15.5% 


78.9% 


6.1% 



Note. — Same as Table but where all 
areas with \b\ < 30° have been included in 
the mask. 



*X 2 computed only for thresholds between 
— 2.5o"o and 2.5 cjq. 



2004), allowing for non-zero shifts and arbitrary ampli- 
tudes. The results from these computations are shown 
in Tables |S] (KpO base mask) and (larger base mask) . 
Only the amplitudes and the asymmetry parameters are 
shown; the shift parameters are in all cases well inside 
the 2cr range, except for the smallest scales (0?53 and 



0?64 FWHM) . The number in parentheses here shows 
the percentage of the simulated maps having the fitted 
parameter smaller than for the WMAP data. 

We see that the genus amplitude A (and especially 
when fitted for v in the range — 2.5tro to — 0.2<7o, i.e., 
A J) is very high in the northern hemisphere and rela- 
tively small in the southern hemisphere, compared with 
the simulations. In particular, the 3?40 scale once again 
stands out as special (2?98 scale with the |6| > 30° base 
mask); in the northern hemisphere the negative thresh- 
old amplitude is so large that only 1 out of 5000 simula- 
tions has a larger amplitude. Using the larger base mask, 
no simulations have as large A_ as the WMAP map at 
2?98 FWHM. In adition, the asymmetry parameter Ag is 
clearly large at the 2 a level in the northern hemisphere. 

In the southern hemisphere, the genus amplitudes are 
relatively, but not extremely, low. However, here the 
asymmetry parameter is rather small, at the 99% level, 
for 2? 98 smoothing and use of the KpO base mask. This 
result is in fact quite peculiar, and it is difficult to decide 
how much significance one should attach to it; it is very 
sensitive to smoothing scale, and it is not seen in the 
larger base mask measurement. One could suspect that 
this result is caused by some features near the Galactic 
plane. However, we have not been able to locate any re- 
gion that by exclusion brings the results to an acceptable 
level. One example from this search is shown in Table 
marked by an asterisk. In this case we have removed two 
disks of 30° radius centered on (I, b) = (330°, -10°) and 
(200°, —20°), correspond ing to the two large cold spots 
discussed bv iParkl i|2004fl . Obviously, the peculiar 99% 
result is not connected to these regions. 

The genus for both the northern and the southern 
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TABLE 5 

Genus results, KpO mask 



FWHM 



Full sky 



0?53 
0?64 
0?85 
1?28 
1?70 
2? 13 
2? 55 
2? 98 
3? 40 
3? 83 
4? 26 



0.087 
0.092 
0.091 
0.082 
0.078 
0.080 
0.082 
0.081 
0.095 
0.096 
0.094 



(87.3%) 
(87.4%) 
(91.5%) 
(91.9%) 
(90.4%) 
(95.3%) 
(96.4%) 
(93.6%) 
(99.5%) 
(99.1%) 
(96.9%) 



0.086 
0.090 
0.089 
0.080 
0.076 
0.080 
0.083 
0.083 
0.101 
0.104 
0.097 



(82.6%) 
(79.3%) 
(86.4%) 
(85.4%) 
(84.2%) 
(93.8%) 
(96.0%) 
(92.9%) 
(99.6%) 
(99.5%) 
(96.1%) 



0.091 
0.096 
0.094 
0.087 
0.084 
0.080 
0.080 
0.083 
0.081 
0.081 
0.091 



(98.0%) 
(98.9%) 
(98.3%) 
(98.4%) 
(98.6%) 
(94.4%) 
(92.2%) 
(93.8%) 
(86.8%) 
(82.0%) 
(93.9%) 



-0.024 
-0.034 
-0.027 
-0.038 
-0.049 
0.000 
0.018 
-0.001 
0.111 
0.128 
0.031 



( 5.5%) 
( 2.3%) 
( 8.1%) 
(10.3%) 
(11.3%) 
(53.9%) 
(67.5%) 
(55.2%) 
(93.5%) 
(93.2%) 
(67.7%) 



Northern hemisphere 



0?85 
1?28 
1?70 
2? 13 
2? 55 
2? 98 
3? 40 
3? 83 
4? 26 



0.100 (97.5%) 
0.090 (93.9%) 
0.091 (97.9%) 
0.095 (98.7%) 
0.104 (99.9%) 
0.110 (99.9%) 
0.122 (99.9%) 
0.116 (99.5%) 
0.113 (97.7%) 



0.099 (96 
0.089 (91 
0.091 (96. 
0.097 (98 
0.110 (99 
0.120 (99 
0.134 (99 
0.129 (99 
0.123 (98 



9%) 
0%) 
6%) 

S) 
9%) 
9%) 
9%) 
7%) 
0%) 



0.104 
0.094 
0.092 
0.087 
0.085 
0.085 
0.092 
0.086 
0.094 



(99.7%) 
(98.4%) 
(97.7%) 
(89.4%) 
(81.0%) 
(76.8%) 
(85.7%) 
(65.8%) 
(77.3%) 



-0.021 
-0.027 
-0.004 
0.056 
0.132 
0.170 
0.188 
0.204 
0.136 



(18.3%) 
(20.4%) 
(47.3%) 
(85.4%) 
(97.6%) 
(98.3%) 
(97.5%) 
(97.0%) 
(87.2%) 



Southern hemisphere 



0?85 
1?28 
1?70 
2? 13 
2? 55 
2? 98 
2? 98* 
3? 40 
3? 83 
4? 26 



0.084 (23.3%) 
0.077 (38.6%) 
0.070 (26.1%) 
0.069 (28.7%) 
0.068 (27.0%) 
0.058 ( 7.6%) 
0.064 (16.0%) 
0.078 (54.2%) 
0.082 (58.6%) 
0.086 (63.2%) 



0.081 
0.074 
0.066 
0.065 
0.062 
0.047 
0.054 
0.082 
0.083 
0.083 



(15.5%) 
(26.6%) 
(15.9%) 
(19.4%) 
(16.8%) 
( 2.4%) 
( 6.5%) 
(63.3%) 
(60.7%) 
(55.4%) 



0.088 
0.081 
0.078 
0.076 
0.077 
0.083 
0.094 
0.067 
0.077 
0.089 



(46.3%) 
(58.9%) 
(60.7%) 
(54.0%) 
(57.9%) 
(69.9%) 
(89.2%) 
(20.3%) 
(43.2%) 
(67.8%) 



-0.039 ( 5.6%) 
-0.044 (10.7%) 
-0.083 ( 4.4%) 
-0.076 (12.4%) 
-0.108 ( 9.2%) 
-0.271 ( 0.4%) 
-0.269 ( 0.9%) 
0.100 (86.0%) 
0.036 (65.4%) 
-0.038 (43.3%) 



Note. — The fitted amplitude of the genus for the full range in v, 
for negative u, and for positive u, and the genus asymmetry parame- 
ter, using the co-added WMAP map with KpO mask. The numbers in 
parentheses show the percentage of the 5000 simulations with a value 
smaller than for the WMAP map. 
*The two large cold spots discussed by Park (2004) are masked out. 



TABLE 6 

Genus results, KpO mask, using only |fe| > 30° 



FWHM 



A;/ 



Northern hemisphere 



2? 13 
2? 55 
2? 98 
3? 40 
3? 83 
4? 26 



0.106 (99.4%) 
0.115 (99.8%) 
0.131 (99.9%) 
0.125 (99.5%) 
0.117 (95.4%) 
0.120 (93.4%) 



0.110 (99.4%) 
0.124 (99.9%) 
0.144 ( 100%) 
0.137 (99.6%) 
0.129 (96.2%) 
0.131 (92.8%) 



0.094 
0.091 
0.095 
0.101 
0.099 
0.106 



(88.8%) 
(76.7%) 
(79.8%) 
(83.0%) 
(73.2%) 
(77.2%) 



0.079 (88.6%) 
0.155 (97.0%) 
0.206 (97.8%) 
0.151 (90.8%) 
0.132 (83.5%) 
0.107 (75.6%) 



Southern hemisphere 



2? 13 
2? 55 
2? 98 
3? 40 
3? 83 
4? 26 



0.078 (43.7%) 
0.079 (46.2%) 
0.068 (15.3%) 
0.077 (33.6%) 
0.108 (89.1%) 
0.104 (79.1%) 



0.077 (40.9%) 
0.079 (44.5%) 
0.064 (14.1%) 
0.082 (45.6%) 
0.127 (96.6%) 
0.123 (90.8%) 



0.084 (64.0%) 
0.081 (48.5%) 
0.077 (32.8%) 
0.065 (10.3%) 
0.064 ( 8.4%) 
0.069 (13.0%) 



-0.044 (27.7%) 
-0.014 (45.7%) 
-0.088 (21.2%) 
0.117 (85.6%) 
0.334 (99.6%) 
0.279 (97.5%) 



Note. — Same as Table QjJ but with all areas with |6| < 30° added 
to the mask. 
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Negative threshold amplitude, A_ 

Fig. 7. — Results from measurements of j4^orth at 30 4Q p-\yHM. 
The best-fit Gaussian distribution has a mean of 0.0771 and a 
standard deviation of 0.0157; therefore, the observed WMAP value 
of 0.134 (thin vertical line) formally corresponds to a 3.6 c effect. 



hemispheres at 2? 98, 3? 40, and 3? 83 FWHMs is shown 
in Figure EJ Here we clearly see the origin of the ex- 
treme numbers found above: the northern hemisphere 
genus at around 3? 40 F WHM lies well outside the lower 
2 cr confidence region and even outside the 3 a region for 
an extended range of thresholds. In addition, the genus 
has a very small amplitude at negative thresholds on the 
southern hemisphere, but an average amplitude at posi- 
tive thresholds, resulting in a rather extreme asymmetry 
parameter. 

The distribution of negative threshold amplitudes A- 
for the northern hemisphere at 3? 40 FWHM is shown 
for the 5000 simulations in Figure together with the 
observed WMAP value of 0.134 (vertical line). The dis- 
tribution is close to Gaussian with a mean of 0.0771 and a 
standard deviation of 0.0157. Thus, the observed WMAP 
A_-value is formally larger than the simulations at 3.6 a. 
However, this number is likely to be slightly overesti- 
mated: since the amplitude is a strictly positive number, 
we should expect the true distribution to be somewhat 
skewed toward large values. On the other hand, this 
effect cannot be very significant in our case, given the 
good fit between the histogram and the Gaussian ap- 
proximation. The case seems solid for rejecting the null 
hypothesis at considerably more than 3 a confidence. 

6.3. Comparison with alternative maps 

So far we have only studied the co-added WMAP 
map, which is dominated by the Q-band. The influ- 
ence of residual foregrounds could therefore be an im- 
portant concern when interpreting these results. In or- 
der to study this issue closer, we now reestimate the 
northern hemisphere statistics (using the KpO mask), 
this time for a set of seven different maps: (1) the 
co-added WMAP map, (2) the Tegmark et al. cleaned 
map (Tegmark, de Oliveira-Costa, & Hamilton 2003), 
3) the WMAP Intern al Linear Combination (ILC) map 
Bennett et aTll2003bJ) . (4)-(6) the averaged Q-, V- and 
VF-bands, and (7) a synchrotron-corrected map [defined 
by [2.65 Ka- KJ/1.65; see Vielva et al. 2004)]. 
Each of these maps has its own effective beam profile, 



and in order for a direct comparison to be meaningful, 
we have to pre-smooth each map so that they have a 
common resolution. This was achieved by first decon- 
volving the old beam and subsequently convolving with 
a 1° Gaussian beam. The only exception is the ILC map, 
which is already smoothed to the appropriate resolution. 

The simulations that we carry out in this paper are 
very CPU-intensive, and it would therefore be very time 
consuming to perform a complete analysis of each of the 
above-mentioned maps. However, after pre-smoothing 
each map to a common resolution, only the effect of noise 
can potentially modify the results, and on the large scales 
that we consider, this effect is very small. We may there- 
fore instead use the existing simulations when interpret- 
ing the new results. 

The two panels of Figure 03 show the results from this 
analysis for the spectral parameter 7 and the negative 
threshold genus amplitude A_ . Each panel contains two 
sets of fundamentally different elements. First, the gray 
1, 2, and 3<r confidence regions and the light green curve 
with green filled circles (the co-added WMAP map) show 
the results of ^6.21 and are thus not pre-smoothcd to 1° 
resolution. The seven remaining curves correspond to 
the maps listed above and are pre-smoothed. 

The effect of the pre-smoothing operation may be re- 
constructed by comparing the two green curves (with cir- 
cles and with diamonds), which show the results for the 
co-added WMAP map, un-smoothed and smoothed, re- 
spectively. For the spectral parameter 7 (left panel), we 
see that the pre-smoothing process effectively suppresses 
7 on small scales, but increases it very slightly on large 
scales. However, this effect is both small and stable on 
large scales, and is thus not difficult to account for. The 
effect on the genus amplitude (right panel) is slightly 
less predictable but follows the same pattern; the pre- 
smoothed function lies generally a little higher than the 
un-smoothed function on large scales. 

The main conclusion to be drawn from these plots, 
however, is that the quantities discussed above show 
a remarkable stability with respect to both frequency 
and foreground correction method. If we neglect the 
synchrotron-corrected map and the WMAP ILC map, 
we see that the scatter in 7 is much less than la, as 
seen by the width of the confidence regions. The same 
is true for the genus amplitude on the scale of most in- 
terest, namely, 2?98 FWHM. We note that the WMAP 
team warns (see the LAMBDA Web site) that the ILC 
map should not be used for CMB studies, and indeed it 
is the most deviant with respect to the other maps in the 
analysis. 

The final issue to consider is the effect of point sources. 
We therefore now apply a median filter to the co-added 
map, compute the same set of quantities for this filtered 
map, and compare the results to those of the unfiltered 
map. Our median filter is implemented by replacing each 
pixel value that is excluded by the WMAP point-source 
mask, i.e., the mask consisting of about 700 disks of 0?5 
radius, by the median evaluated over a 1° disk around 
that pixel. Such filters are commonly used in, e.g., radio 
astronomy, to eliminate outliers. 

The results from this experiment are shown in Fig- 
ure El Once again, we plot 7 and A- estimated in the 
northern hemisphere (KpO mask), since these quantities 
represent the two strongest detections presented in this 
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Fig. 8. — Scale dependence of the spectral parameter 7 (left) and the negative threshold genus amplitude A— (right), computed on the 
northern hemisphere for seven different maps. The gray bands show the 1, 2, and 3cr confidence regions computed from simulations. All 
maps have been pre-smoothed to a common 1° resolution, except for the co-added WMAP map, which is plotted both with and without 
pre-smoothing. The confidence bands are computed from simulations that are not pre-smoothed. 
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Fig. 9. — Effect of point sources on the spectral parameter 7 (left) and the negative threshold genus amplitude A— (right). The functions 
are computed on the northern hemisphere, restricted to the KpO base mask. Circles indicate the functions computed from the raw maps, 
and squares indicate the functions computed from median-filtered maps. 



paper. The results are shown for two different maps, the 
raw co-added map and the WMAP template-corrected 
map. In both cases the effect of the median filter may 
be seen by comparing the solid and dashed lines (circles 
and squares), which correspond to the raw map and the 
median filtered map, respectively. 

It is apparent that median filtering changes 7 by a 
negligible amount on large scales while on small scales 7 
is very slightly decreased. The effect on the genus am- 
plitude A- is less predictable, as the two curves lie both 
above and below each other. However, it is interesting to 
note that the filtered maps are associated with a slightly 
larger amplitude than the raw maps at 2? 98 FWHM and 
a slightly smaller amplitude at 3?83 FWHM. Of course, 
it is difficult to estimate the significance of this. In any 
case, it is clear that the impact of point sources is too 
small to explain the reported large-scale results in this 
paper. 



6.4. Correlations between the measurements 

In the previous sections we presented measurements of 
several important statistics. In this section we go one 
step further and seek to understand the correlations be- 
tween some of these statistics. This goal is important 
both in order to understand the nature of the quantities 
themselves and for assessing the combined significance of 
the quoted results. In the following we focus on the full- 
sky measurements smoothed with a f?28 FWHM Gaus- 
sian beam. 

Figure ^| shows the following scalar quantities plotted 
pair-wise against each other: the temperature standard 
deviation cr versus the spectral parameter 7, the spectral 
parameter 7 versus the full range genus amplitude A, 
and the negative threshold genus amplitude A_ versus 
the positive threshold amplitude A + . The values from 
each simulated realization are marked by a small filled 
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Fig. 10. — Correlation plots of scalar quantity pairs, computed from maps smoothed with a 1?28 FWHM beam and constrained to the 
KpO base mask. Each small filled circle corresponds to 1 of 2000 Monte Carlo realizations, and the cross marks the WMAP value. The 
plotted relationships are (a) the temperature standard deviation cfq versus the spectral parameter 7, (fe) the spectral parameter 7 versus 
the full range genus amplitude A, and (c) the negative threshold genus amplitude A— versus the positive threshold amplitude A+. 



circle, and the WMAP results are marked by a cross. 

In Figure ITUb we see that the correlation between <to 
and 7 is very strong and negative. Of course, given that 
7 = a\j '(<7n(72), this is no surprise. Nevertheless, it is 
well worth noticing that the spectral parameter 7 may in 
many respects be identified with the standard deviation 
of the map. 

Figure ITUb demonstrates the clear correlation between 
the spectral parameter 7 and the full range genus am- 
plitude A. Although the scatter is non-negligible, a high 
7-value is probably accompanied by a large genus am- 
plitude. Thus, given the large 7-value of WMAP, the 
large observed genus amplitude seen in both Table 0and 
Figure [S] is not unexpected. 

Finally, in Figure 110b we see that there is only a very 
weak correlation between the negative threshold genus 
amplitude A_ and the positive threshold amplitude A + . 

The conclusions to be drawn from these plots are there- 
fore the following: 

• The spectral parameter 7 is almost inversely pro- 
portional to the standard deviation of the map. 

• The spectral parameter 7 and the genus amplitudes 
are correlated, and the corresponding numbers in 
Tables and Table may not be considered as 
independent results. 

• The negative and positive threshold genus ampli- 
tudes are only very mildly correlated. 

Given these relations, we may now make a few con- 
nections between the results presented in H6.2I First, 
we noted that the WMAP spectral parameter 7 was 
large in both the full-sky measurements and particularly 
the northern hemisphere. However, it was not signif- 
icantly different from the simulations in the southern 
hemisphere. These results are thus in perfect agreement 
with the genus measurements shown in Figures and |H| 
and the parameters given in Table — both 7 and the 
genus amplitudes are large in the northern hemisphere, 
but fully acceptable in the southern hemisphere. The 



measurements of 7 and the genus amplitudes may obvi- 
ously not be regarded as independent. 

One further interesting connection may be made 
through the relation between 7 and the standard devia- 
tion en, which may be taken as a very crude measure of 
the overall power in the map, determined primarily by 
the large-scale modes. Since 7 is very large in the north- 
ern hemisphere, the standard deviation is very small; in 
other words, there is little power in this re gion. This re- 
sult th erefore supports the conclusions of lEriksen et alJ 
( 2004) , that there is a clear lack of large-scale power in 
the northern hemisphere. Conversely, the genus results 
are generally stronger than the 7 results, and given the 
relatively high Ag results, a complete explanation must 
probably be formulated as a combination of power spec- 
trum and other non-Gaussianity effects. 

7. CONCLUSIONS 

In this paper we have presented the three two- 
dimensional Minkowski functionals and the length of the 
skeleton of the co- added WMAP map, as well as a num- 
ber of other statistics, and compared these with sim- 
ulations based on Gaussian perturbations with a best- 
fit power spectrum and WMAP-specitic noise and beam 
properties. On scales between 1° and 4°, all the statis- 
tics except for the genus Minkowski functional accept the 
model. On scales smaller than ~ 0?9 there are too many 
stationary points, a result that is interpreted as an effect 
of point sources. 

By studying each Galactic hemisphere separately, it 
was found that on scales around 3?4, the northern hemi- 
sphere deviates significantly from the simulations. Most 
importantly, the genus functional has a very large ampli- 
tude on these scales; only 1 Monte Carlo realization out 
of 5000 simulations has a larger amplitude for negative 
thresholds. Interestingly, the southern hemisphere genus 
has a fairly small amplitude around the same scales. 
Furthermore, the spectral parameter 7 is higher in the 
WMAP data than in the simulations on scales ~ 3° in the 
northern hemisphere. These asymmetries may be com- 
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pared wit h thos e rece ntly announced by lEriksen et alJ 
( 2004) and lParkl <|2004fl . Park also reports an asymmetry 
in the genus amplitudes between the northern and south- 
ern hemispheres. However, these results are found on sig- 
nificantly smaller scales than those considered here, and 
it is therefore difficult to estab l ish a direct link between 
our findings and those of lParkl l)2004(l . 

Many of the results presented in this paper depend on 
both the power spectrum and the higher-order statistical 
properties of the field. For example, the genus ampli- 
tude depends to a large extent on the power spectrum, 
and the strong results found on the northern hemisphere 
may therefore be influenced by both power spectrum and 
non-Gaussianity issues. However, the asymmetry param- 
eter, Ag, is a much cleaner test of non-Gaussianity, and 
we found that this was extreme at the 99% level in the 
southern hemisphere in a few cases. This result may 
therefore be t aken as support for the results presented by 
iVielva et alJ |200l . who find non-Gaussian signatures in 
the southern hemisphere using spherical wavelets. How- 
ever, the WMAP Ag parameter is also large at the 2 a 
level in the north aroun d 3° FWHM, and th is is not sup- 
ported by the results of IVielva et al.1 l)2004[L It is there- 
fore quite difficult to determine if our results support 
the previously reported non-Gaussian effects, or whether 
these findings are independent. 

It is in any case clear that the currently accepted Gaus- 
sian model has problems accounting for the statistical 



properties of the WMAP data on large and intermediate 
scales. In this paper we have focused particularly on the 
co-added Q + V + W WMAP data, but have also shown 
that the results are very stable with respect to frequen- 
cies and different foreground correction methods. Thus, 
uncertainties in foreground contributions or residuals are 
unlikely to resolve these issues. While it certainly is too 
early to conclude that these detections are of primor- 
dial origin, the possibility should certainly be taken into 
account. The 2-year WMAP data, and eventually the 
Planck data, will be essential in determining the nature 
of these detections. 
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